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ABSTRACT 

Large-scale polarization of the cosmic microwave background measured by the WMAP 
satellite requires a mean optical depth to Thomson scattering, r e ~ 0.17. The reion- 
ization of the universe must therefore have begun at relatively high redshift. We have 
studied the reionization process using supercomputer simulations of a large and repre- 
sentative region of a universe which has cosmological parameters consistent with the 
WMAP results. Our simulations follow both the radiative transfer of ionizing photons 
and the formation and evolution of the galaxy population which produces them. A 
previously published model with a standard stellar IMF, with a moderate photon es- 
cape fraction from galaxies and with ionizing photon production as expected for zero 
metallicity stars produces r e = 0.104, which is within 1.0 to 1.5a of the "best" WMAP 
value. Values up to 0.16 can be produced by taking larger escape fractions or a top 
heavy IMF. The data do not require a separate populations of "miniquasars" or of 
stars forming in objects with total masses below 10 9 M Q . Reconciling such early reion- 
ization with the observed Gunn-Peterson troughs in z > 6 quasars may be challenging. 
Possible resolutions of this problem are discussed. 

Key words: galaxies: formation - cosmology: theory - large-scale structure - cosmol- 
ogy: observations 



MOTIVATION 



It has recently become possible to study the reionization of 
the universe in some detail. Massive numerical computations 
can now follow not only the clustering of dark matter and 
gas, but also the formation of galaxies and the propagation 
of ionizing photons in a highly inhomogeneous, partially ion- 
ized, intergalactic medium (IGM). In spite of many poorly 
understood details concerning the physics of star forma- 
tion, and the approximations inherent in the various numer- 
ical treatments of radiative transfer (for a recent review see 
Maselli, Ferrara & Ciardi 2003), a number of independent 
studies have converged on a relatively late (z r ^ 8 to 10) 
epoch for complete reionization of the IGM within current 
"concordance" (i.e. flat, A-dominated) cosmological models 
(e.g. Ciardi et al. 2000, hereafter CFGJ; Gnedin 2000; Ra- 
zoumov et al. 2002; Ciardi, Stoehr & White 2003, hereafter 
CSW). 

This conclusion appears challenged by results from the 
WMAP satellite (Kogut et al. 2003; Spergel et al. 2003). 
This experiment has detected an excess in the CMB TE 
cross-power spectrum on large angular scales [t < 7) indi- 
cating an optical depth to the CMB last scattering surface 
of r e = 0.16. The uncertainty quoted for this number de- 
pends on the analysis technique employed. Fitting the TE 



cross power spectrum to ACDM models in which all param- 
eters except r e take their best fit values based on the TT 
power spectrum, Kogut et al. (2003) obtain a 68% confi- 
dence range, 0.13 < r e < 0.21. Fitting all parameters simul- 
taneously to the TT+TE data, Spergel et al. (2003) obtain 
0.095 < r e < 0.24. Including additional data external to 
WMAP, these authors were able to shrink their confidence 
interval to 0.11 < r e < 0.23. Finally, by assuming that the 
observed TT power spectrum is scattered to produce the 
observed TE cross-power spectrum Kogut et al. (2003) infer 
0.12 < r e < 0.20. It is unclear which of these uncertainty 
estimates to prefer. Most r e values in these ranges require 
a substantial fraction of the universe to be ionized before 
redshift 10. For example, 0.12 < r e < 0.20 translates into 
13 < z r < 19 in a model where reionization is instantaneous. 

Apparently, reionization occurred earlier than expected. 
Is this discrepancy real? In practice, recovering the "reion- 
ization epoch" is subject to some ambiguity related to the 
specific reionization history assumed (see e.g. Bruscoli, Scan- 
napieco & Ferrara 2002). What might reionization models 
have overlooked? The discrepancy is not dependent on the 
particular cosmological parameters adopted, since WMAP 
has confirmed the previous concordance model. Hence the 
"oversight" must be of astrophysical nature. Several effects 
might have produced rapid early evolution: a contribution 
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from low-mass pregalactic systems (in the jargon, Popu- 
lation III objects); unexpectedly high star formation effi- 
ciency; unexpectedly high ionizing photon production; an 
unexpectedly large probability for ionizing photons to escape 
into the IGM; a possible population of early "miniquasars" . 
In this paper, we point out that efficient production and es- 
cape of ionizing photons is sufficient to account for the data 
within conventional galaxy formation models. The data do 
not require very massive stars or miniquasars, although the 
high efficiencies needed may point to near-zero metallicities 
or to a top-heavy stellar Initial Mass Function (IMF) at 
early times. 



2 SIMULATING COSMIC REIONIZATION 

In CSW we studied cosmic reionization through a combina- 
tion of high-resolution N-body simulations (to describe the 
distribution of dark matter and diffuse gas), a semi-analytic 
model of galaxy formation (to track the sources of ioniza- 
tion) and the Monte Carlo radiative transfer code CRASH 
(to follow the propagation of ionizing photons in the IGM; 
Ciardi et al. 2001, CFMR; Maselli, Ferrara & Ciardi 2003). 
Here, we briefly summarize the features of the above ap- 
proach which are relevant to the present study. 

The simulations are based on a ACDM "concordance" 
cosmology with Q m —0.3, CIa~0-7, h =0.7, fi(,=0.04, n—1 
and as— 0.9. These parameters are within the WMAP exper- 
imental error bars (Spergel et al. 2003). The re-simulation 
technique described in Springel et al. (2001, hereafter 
SWTK) and the N-body code GADGET (Springel, Yoshida 
& White 2001) were used to follow at high resolution the 
dark matter distribution within an approximately spheri- 
cal subregion of diameter about 50ft" 1 Mpc within a much 
larger cosmological volume. Within this subregion a cube of 
comoving side L = 20ft" 1 Mpc was used for the detailed 
radiative transfer modeling. The location and mass of dark 
matter halos was determined with a friends-of-friends algo- 
rithm. Gravitationally bound substructures were identified 
within the halos with the algorithm SUBFIND (SWTK) and 
were used to build the merging tree for halos and subhalos 
following the prescription of SWTK. The smallest resolved 
halos have masses of M ~ 10 9 M and they start forming in 
our box at z ~ 20. We model the galaxy population via the 
semi-analytic technique of Kauffmann et al. (1999) as imple- 
mented by SWTK. At the end of this process we obtain a 
mock catalogue of galaxies for each of the 51 simulation out- 
puts, containing for each galaxy, among other quantities, its 
position, stellar mass and star formation rate. The simula- 
tions match the observed global star formation rate density 
evolution for z < 5. The N-body simulation used in this pa- 
per is the M3 simulation of CSW and all galaxy formation 
modeling is identical to that in the earlier paper. 

Reionization simulations require a mass resolution high 
enough to follow the formation and evolution of the objects 
producing the bulk of the ionizing radiation. At the same 
time, large simulation volumes must be considered to avoid 
biases due to cosmic variance on small scales. Moreover, a re- 
liable treatment of the radiative transfer of ionizing photons 
in the IGM is needed. So far, although several numerical ap- 
proaches have been proposed (e.g. Gnedin & Ostriker 1997; 
CFGJ; Chiu & Ostriker 2000; Gnedin 2000; Razoumov et 



RUN 


IMF 


fesc 


S5 


Salpctcr 


5% 


S20 


Salpctcr 


20% 


L20 


Larson 


20% 



Table 1. Parameters of the simulations: Initial Mass Function, 
IMF; photon escape fraction, fesc- 



al. 2002) only the simulations described in CSW fulfill all 
the above requirements (refer to that paper for an extensive 
discussion) . 

One of the aims of the present study is to assess the 
influence of the IMF on reionization. We have thus inferred 
the emission properties of our model galaxies by assuming 
a time-dependent spectrum of a simple stellar population of 
metal-free stars in the mass range up to 40 M Q , with two 
different IMFs (CFMR): (i) a Salpeter IMF and (ii) a Larson 
IMF, i.e. a Salpeter function at the upper mass end which 
falls off exponentially below a characteristic stellar mass, 
M c — 5 Mq. The Larson IMF has a specific photon flux at 
the Lyman continuum which is about 4 times the Salpeter 
one. The total number of ionizing photons per solar mass 
for a Larson IMF is « 4 x 10 60 M" 1 . 

Of the emitted ionizing photons, only a fraction / esc 
will actually be able to escape into the IGM. This quantity 
is poorly determined both theoretically and observationally: 
actually, f eac may well vary with, e.g., redshift, mass and 
structure of a galaxy, as well as with the ionizing photon 
production rate (e.g. Wood & Loeb 1999; Ricotti & Shull 
2000; Ciardi, Bianchi & Ferrara 2002). However, at moderate 
redshift, recent results on the opacity evolution of the Lya 
forest (Bianchi, Cristiani & Kim 2001) constrain f e3C to be 
smaller than 20% in order not to over-produce the cosmic 
UV background. Although this limit may not apply to the 
high redshift universe, we consider this value a physically 
plausible upper limit. 

Given these assumptions, the simulations of reionization 
described above have been run for the three different param- 
eter combinations in Table 1. Run S5 (extensively described 
in CSW as the M3 case), yields the lowest IGM ionization 
power input, whereas L20 maximizes it; run S20 is an in- 
termediate case. The critical parameter differentiating these 
runs is the number of ionizing photons escaping a galaxy 
into the IGM for each solar mass of long-lived stars which it 
forms. Each simulation can be thought of as corresponding 
to models with varying IMF, stellar luminosity beyond the 
Lyman limit and value of f esc , provided this quantity is kept 
fixed (see CSW). 

We quantify the agreement between our simulations and 
the WMAP data primarily through the optical depth to elec- 
tron scattering, r e , given by: 



Jo 



a T n e (z')c 



dt 



dz' 



dz' , 



(1) 



where <jt = 6.65 x 10~ 25 cm 2 is the Thomson cross sec- 
tion and rie(z') is the mean electron number density at z' . 
We will compare our simulated r e with the "model indepen- 
dent" estimate of Kogut et al. (2003) r e = 0.16 ± 0.04 (68% 
confidence range), but it should be borne in mind that the 
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Figure 1. Slices through the simulation boxes. The six panels show the neutral hydrogen number density for the L20 (upper panels) 
and the S20 (lower panels) runs, at redshifts, from left to right, z = 17.6, 15.5 and 13.7. The box for the radiative transfer simulation has 
a comoving length of L = 20/i -1 Mpc. 



Spergel et al. (2003) joint analysis of WMAP and external 
data suggests a significantly larger allowed range. 



3 RESULTS 

In order to clarify the differences between our runs, we start 
from a visual inspection of simulated maps. Fig. 1 shows 
the redshift evolution of the Hi number density for the L20 
(upper panels) and the S20 (lower) runs (illustrative maps 
for the S5 runs can be found in CSW). Highly ionized re- 
gions (dark areas) are produced by the young galaxies in 
the box and are well resolved in the maps. They initially oc- 
cupy a small fraction of the volume, and are typically larger 
in the L20 run because of the higher ionizing power of the 
sources. Their shape, particularly for larger ones, appears 
distorted by nearby high density peaks. (To a good approx- 
imation these correspond to peaks in the Hi distribution.) 
The ionization front slows when it encounter such overdensi- 
ties because of their higher recombination rate. By redshift 
z = 15.5 several bubbles are close to overlap in the L20 run 
(see upper-right corner of the central panel) whereas in the 
S20 run the filling factor is still small. Finally, by z — 13.7 
the overlapping fronts have cleared out most of the volume 
in L20 with tiny Hi islands surviving thanks to their high 
density; reionization in the S20 run, on the other hand, is 
far from complete. 

These remarks can be complemented by the more quan- 
titative analysis shown in Fig. 2, which plots the (normal- 
ized) distribution function of the simulated electron density 
at various redshifts for the two extreme runs, L20 and S5. 



This distribution is defined as the fraction of the total vol- 
ume filled with gas with free electron density in each loga- 
rithmic bin of rie- Its bimodal shape reflects the two-phase 
(neutral + ionized) structure of the IGM. The mean parti- 
cle density in the redshift range 13.3 < z < 18.5 is around 
10~ 3 cm~ 3 , so the rightmost peak has n e ~ n and corre- 
sponds to the ionized phase. Conversely, the left peak is as- 
sociated with mostly neutral gas. The general evolutionary 
trend is a continuous transfer of matter from the left peak to 
the right one as redshift decreases. By z = 13.3 we find that 
the mass-averaged neutral hydrogen fraction is only 2% for 
the L20 model, whereas for the S5 and S20 cases it is still 
89% and 63% respectively. 

The redshift evolution of the volume-averaged ioniza- 
tion fraction, x v , essentially coincident with mass-averaged 
curve (CSW), is reported for each of our three runs in the 
inset of Fig. 3. An L20 vs S20 comparison shows that the 
former typically has an x v value fa 3 times higher, reaching 
complete ionization (x v ~ 1) at z r » 13. Run S20, on the 
other hand, reaches complete reionization at z r fv 11. In run 
S5 reionization proceeds much more gradually and it is only 
at z T ~ 8 that the value x v ~ 1 is reached. 

Finally, we have calculated the evolution of r e (Fig. 3) 
corresponding to the above reionization histories as follows. 
Prior to complete reionization, n e (z) in eq. 1 is obtained 
from the simulations; after the reionization epoch, we sim- 
ply assume complete H and Hei ionization throughout the 
box. We also assume Hen reionization at z = 3. The three 
runs yield the values r e = 0.104 (S5), r e = 0.132 (S20) 
and T e = 0.161 (L20). A value r e = 0.16 is also obtained if 
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10- 7 io- 6 10- 5 10- 4 10- 3 io- 2 

n e [cm" 3 ] 



Figure 2. Evolution of the (normalized) electron density distri- 
bution function, n e , for the L20 (upper panel) and S5 (lower) 
runs. Curves refer to different redshifts, as in the label. 

one assumes instantaneous reionization at z r ~ 16 (dotted 
line), i.e. three redshift units higher than the actual epoch of 
complete reionization in the model. At z — 16 the ionization 
fraction is x v ~ 0.3 in L20. 



4 DISCUSSION 

Our supercomputer simulations of galaxy formation and of 
the propagation of ionization fronts have shown that recent 
WMAP measurements of the optical depth to electron scat- 
tering (r e =0.16±0.04 according to the "model independent" 
analysis of Kogut et al. (2003) or r e =0.17±0.06 according 
to the parameter fits of Spergel et al. (2003)) are easily re- 
produced by a model in which reionization is caused by the 
first stars in galaxies with total masses of a few xlO 9 M0. 
In order to get sufficient early ionization, this phase of star 
formation must supply a relatively high number of ionizing 
photons to the IGM for each solar mass of long-lived stars 
which is formed. This requires some combination of a high 
photon escape fraction, a top-heavy IMF, and a high stel- 
lar production rate for ionizing photons, similar, perhaps, 
to that typically inferred for metal-free stars. Among the 
three models we have explored, the "best" WMAP value for 
t,. is matched assuming a moderately top-heavy IMF and 
an escape fraction of 20%. A Salpeter IMF with the same 
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Figure 3. Redshift evolution of the electron optical depth, r e , 
for the S5 (long-dashed line), S20 (short-dashed) and L20 (solid) 
runs. The dotted line refers to sudden reionization at z = 16. 
The shaded region indicates the optical depth r e = 0.16 ± 0.04 
(68% CL) implied by the Kogut et al. (2003) "model independent" 
analysis. In the inset the redshift evolution of the volume-averaged 
ionization fraction, x v , is shown for the three runs. 

escape fraction gives r e = 0.132, which is still within all 
the suggested 68% confidence ranges. Decreasing f esc to 5% 
gives r e = 0.104, which disagrees with WMAP only at the 
1.0 to 1.5a level depending on which uncertainty estimate 
one adopts. All these models assume that ionizing photons 
are produced at the level expected for metal-free stars. It is 
clearly possible to reproduce the experimental data without 
invoking exotica such as very massive stars, early "mini- 
quasars", or minihalos, i.e. halos cooled by molecular hy- 
drogen with virial temperatures < 10 4 K. 

In our best-fitting model (L20) reionization is essen- 
tially complete by z r w 13. This is difficult to reconcile with 
observations of the Gunn-Peterson effect in z > 6 quasars 
(Becker et al. 2001; Fan et al. 2002). These imply a volume- 
averaged neutral fraction above 10 -3 and a mass-averaged 
neutral fraction ~ 1% at z = 6. These values are reached 
in our L20 run at z « 13. Even for our S5 model, they are 
reached before z ~ 8. Thus the neutral fraction at late times 
must be higher than in our models if the Gunn-Peterson data 
are to be consistent with the WMAP findings. A fascinat- 
ing (although speculative) possibility is that the universe 
was reionized twice (Cen 2002; Wyithe & Loeb 2002) with a 
relatively short redshift interval in which the IGM became 
neutral again. The maximum thickness of the neutral layer 
(ending at z ~ 6) for which the L20 run remains consis- 
tent with the WMAP data, is Az ~ 3; this sets the end of 
the first reionization epoch at z w 9. Distinguishing single- 
from two-epoched reionization models is in principle possible 
from the TE cross-correlation spectra (Naselsky & Chiang 
2003); however, because of the high sensitivity required, this 
measurement will have to await the Planck mission. 

What mechanism could have reduced the production 
of ionizing photons enough to allow recombination at the 
end of the first reionization epoch? Suppression of the 
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galaxy formation process itself seems implausible. An in- 
crease of the relevant filtering scale could have suppressed 
galaxy formation in halos below a typical circular velocity 
of 30 - 40 km s" 1 (Gncdin 2000) but this is well below the 
scale of the halos which host most of the star formation in 
our models at redshifts below 10. Feedback by SNe might 
reduce formation efficiencies (e.g. Madau, Ferrara & Rees 
2001; Scannapieco et al. 2000) but seems likely to be more 
effective at the high redshifts where we need efficient star 
formation than at the later times when recombination is 
supposed to occur. A more plausible possibility may be that 
galaxy formation continues apace, but that the efficiency of 
ionizing photon production drops dramatically, perhaps as 
a result of increasing stellar metallicity, of decreasing escape 
fraction, or of changing IMF. 

In many ways our models are quite conservative. We 
have simply extrapolated conventional models for galaxy for- 
mation to higher redshift assuming that feedback and star 
formation efficiencies can be scaled down to systems with 
mass a few xlO 9 M0 and circular velocities ~ 60 km s _1 
at 2 w 15. We have adopted optimistic but not implausible 
values for ionizing photon production efficiency and escape 
fraction. We do not require any stars to form in low mass 
objects or through molecular cooling processes, nor do we 
invoke any non-stellar sources of ionizing radiation. 

Processes we have ignored may nonetheless play a sig- 
nificant role. We do not attempt to model the possible ef- 
fects of (unresolved) minihalos. Such objects rely on hydro- 
gen molecules for cooling, and it is unclear whether they 
can form stars at all. As pointed out by CFGJ, and more 
recently by Wyithe & Loeb (2003), their contribution to 
the ionizing photon budget is in any case expected to be 
negligible. Minihalos could also increase the IGM clumping 
factor, thus enhancing its recombination efficiency and pro- 
viding a potentially important sink for photons. If this ef- 
fect were strong, more extreme assumptions about the IMF 
and the escape fractions would be required to ensure early 
reionization. One might then be forced to invoke very mas- 
sive stars. For example, a 200 M star, ending its life as a 
pair-instability supernova (SN 77 ), produces approximately 
3 x 10 3 '(Z I 'Zq) ionizing photons/H-atom. Stars of this mass 
may cease to form once the metallicity of the parent gas 
cloud exceeds Z « 10~ 5±1 Z Q (Bromm et al. 2001; Schnei- 
der et al. 2002). They would then disappear after they have 
contributed about 0.003 - 0.3 photons/H-atom, hence prior 
to reionization. In addition, one could posit stars outside 
the SN 77 mass range, which end up in very massive black 
holes. This may lead to a star formation conundrum (Schnei- 
der et al. 2002): lacking the necessary heavy element pol- 
lution, only extremely heavy stars would form indefinitely. 
Although a mixture of the two populations, dominated by 
BH-progenitor stars might be favored by observations of 
the near infrared background (Salvaterra & Ferrara 2002; 
Magliocchetti, Salvaterra & Ferrara 2003), our current re- 
sults suggest that the WMAP measurement of r e does not 
require any such exotic objects. 
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